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Abstract 

The interaction between visco-elasto-plastic and adhesive particles is the sub¬ 
ject of this study, where “meso-particles” are introduced, i.e., simplified particles, 
whose contact mechanics is not taken into account in all details. A few examples 
of meso-particles include agglomerates or groups of primary particles, or inhomo¬ 
geneous particles with micro-structures of the scale of the contact deformation, 
such as core-shell materials. 

A simple, flexible contact model for meso-particles is proposed, which allows 
to model the bulk behavior of assemblies of many particles in both rapid and slow, 
quasi-static flow situations. An attempt is made to categorize existing contact mod¬ 
els for the normal force, discuss all the essential mechanical ingredients that must 
enter the model (qualitatively) and finally solve it analytically. 

The model combines a short-ranged, non-contact part (resembling either dry 
or wet materials) with an elaborate, visco-elasto-plastic and adhesive contact law. 
Using energy conservation arguments, an analytical expression for the coefficient 
of restitution is derived in terms of the impact velocity (for pair interactions or, 
equivalently, without loss of generality, for quasi-static situations in terms of the 
maximum overlap or confining stress). 

Adhesive particles (or meso-particles) stick to each other at very low impact 
velocity, while they rebound less dissipatively with increasing velocity, in agree¬ 
ment with previous studies. For even higher impact velocities an interesting sec¬ 
ond sticking and rebound regime is reported. The low velocity sticking is due to 
non-contact adhesive forces, the first rebound regime is due to stronger elastic and 
kinetic energies with little dissipation, while the high velocity sticking is generated 
by the non-linearly increasing, history dependent plastic dissipation and adhesive 
contact force. As the model allows also for a stiff, more elastic core material, this 
causes the second rebound regime at even higher velocities. 

Keywords: Meso-scale particles and contact models, Particle collisions, 
Plastic loading-unloading cycles, Sticking, Adhesive contacts, Cohesive pow¬ 
ders, Elasto-plastic material, Core-shell particles 
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Nomenclature 


mi : mass of i th particle. 

at : Radius of i ,h particle. 

nt r : Reduced mass of a pair of particles. 

5 : Contact overlap between particles. 

Vi : Relative velocity before collision. 

Vf : Relative velocity after collision. 

Vi°° : Relative velocity before collision at infinite separation. 

v f °° : Relative velocity after collision at infinite separation, 

v" : Normal component of relative velocity. 

e : Coefficient of restitution. 

e n : Normal coefficient of restitution. 

Ei : Pull-in coefficient of restitution. 

e 0 : Pull-off coefficient of restitution. 

k : Spring stiffness. 

k\ : Slope of loading plastic branch. 

: Slope of unloading and re-loading elastic branch. 
k c : Slope of irreversible, tensile adhesive branch. 

k p : Slope of unloading and re-loading limit branch; end of plastic regime. 

v p : Relative velocity before collision for which the limit case is reached. 

(j)f : Dimensionless plasticity depth. 

<5 m ax : Maximum overlap between particles during a collision. 

<5max : Maximum overlap between particles for the limit case. 

So : Force free overlap = plastic contact deformation. 

5mi n : Overlap between particles at the maximum (negative) attractive force. 

S c '■ Kinetic energy free overlap between particles. 

Wdiss : Amount of energy dissipated during collision. 

77 : Dimensionless plasticity of the contact. 

p : Adhesivity: dimensionless adhesive strength of the contact. 

X : Scaled initial velocity relative to v p . 

f a : Non-contact adhesive force at zero overlap. 

<5 a : Non-contact separation between particles at which attractive force becomes active. 

IS' : Strength of non-contact adhesive force. 
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1 Introduction 


Granular materials and powders are ubiquitous in industry and nature. For this rea¬ 
son, the past decades have witnessed a strong interest in research aiming for better 
understanding and predicting their behavior in all regimes from flow to static as well as 
the transitions between these states. Especially, the impact of fine particles with other 
particles or surfaces is of fundamental importance. The interaction force between two 
particles is a combination of elasto-plastic deformation, viscous dissipation, and ad¬ 
hesion - due to both mechanical contact- and long ranged non-contact forces. Pair 
interactions that can be used in bulk simulations with many particles and multiple con¬ 
tacts per particle are the focus here, and we use the special, elementary case of pair 
interactions to understand them analytically. 

Different regimes can be observed for collisions between two particles: For ex¬ 
ample, a particle can either stick to another particle/surface or it rebounds, depending 
upon the relative strength of adhesion and impact velocity, size and material various 
material parameters JTj. This problem needs to be well understood, as it forms the ba¬ 
sis for understanding rather complex, many-particle flows in realistic systems, related 
to e.g. astrophysics (dust agglomeration, Saturn’s rings, planet formation) or industrial 
processes (handling of fine powders, granulation, filling and discharging of silos). Par¬ 
ticularly interesting are the interaction mechanisms for adhesive materials such as as¬ 
phalt, ice particles or clusters/agglomerates of fine powders (often made of even smaller 
primary particles). Some of these materials can be physically visualized as having a 
plastic outer shell with a stronger and more elastic inner core. Understanding this can 
then be applied to particle-surface collisions in kinetic spraying, where the solid micro¬ 
sized powder particle is accelerated towards a substrate. In cold spray, bonding occurs 
when impact velocities of particles exceed a critical value, which depends on vari¬ 
ous material parameters IBS- However, for even higher velocities particles rebound 
from the surface 00. Due to the inhomogeneity of most realistic materials, their 
non-sphericity and their surface irregularity, one can not include all these details - but 
rather has to focus on the essential phenomena and ingredients, finding a compromise 
between simplicity and realistic contact mechanics. 

1.1 Contact Models Review 

Computer simulations have turned out to be a powerful tool to investigate the physics 
of particulate systems, especially valuable as experimental difficulties are considerable 
and since there is no generally accepted theory of granular flows. A very popular sim¬ 
ulation scheme is an adaptation of the classical Molecular Dynamics technique called 
Discrete Element Method (DEM) (for details see Refs. IMU)- It involves integrating 
Newton’s equations of motion for a system of “soft”, deformable grains, starting from 
a given initial configuration. DEM can be successfully applied to adhesive particles, if 
a proper force-overlap law (contact model) is used. 

The JKR model Ifl6ll is a widely accepted contact model for adhesive elastic spheres 
and gives an expression for the normal force in terms of the normal deformation. Der- 
jaguin et al. cm suggested that the attractive forces act only just outside the contact 
zone, where surface separation is small, and is referred to as DMT model. An in- 
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teresting approach for dry adhesive particles was proposed by Molerus mm, who 
explained consolidation and non-rapid flow of adhesive particles in terms of adhe¬ 
sive forces at particle contacts. Thornton and Yin |20l compared the results of elastic 
spheres with and without adhesion, a work that was later extended to adhesive elasto- 
plastic spheres ED. Molerus’s model was further developed by Tomas, who intro¬ 
duced a complex contact model l22t[24l by coupling elasto-plastic contact behavior 
with non-linear adhesion and hysteresis involving dissipation and a history (compres¬ 
sion) dependent adhesive force. The contact model subsequently proposed by Lud- 
ing I15II25I works in the same spirit as that of Tomas |23l . only reducing complexity 
by using piece-wise linear branches in an otherwise non-linear contact model in spirit 
(as explained later in this study). In the original version fl5l . a short-ranged force be¬ 
yond contact was mentioned, but not specified, which is one of the issues tackled in 
the present study. Contact details, such as a possible non-linear Hertzian law for small 
deformation, and non-linear loading-unloading hysteresis are over-simplified in Lud- 
ings model, as compared to the model proposed by Tomas |23| . This is partly due to 
the lack of the experimental reference data or theories, but also motivated by the wish 
to keep the model as simple as possible. The model consists of several basic mecha¬ 
nisms, i.e., non-linear elasticity, plasticity and adhesion as relevant for, e.g. core-shell 
materials or agglomerates of fine, dry primary powder particles Il26l[27l l. A possible 
connection between the microscopic contact model and the macroscopic, continuum 
description for adhesive particles was recently proposed by Luding et al. l28l . as fur¬ 
ther explored by Singh et al. H291I301I for dry adhesion, by studying the force anisotropy 
and force distributions in steady state bulk shear in the ID, which is further generalized 
to wet adhesion by Roy et al. 11331 1 . or studied under shear-reversal 13411351 . 

Jiang et al. [36 ] experimentally investigated the force-displacement behavior of 
idealized bonded granules. This was later used to study the mechanical behavior of 
loose cemented granular materials using DEM simulations lf37t . Kempton et al. 1381 
proposed a meso-scale contact model combining linear hysteretic, simplified JKR and 
linear bonding force models, to simulate agglomerates of sub-particles. The phe¬ 
nomenology of such particles is nicely described by Dominik and Tielens 1268 . Walton 
et al. H39II40I also proposed contact models in similar spirit as that of Luding |fl5l and 
Tomas 1238 . separating the pull-off force from the slope of the tensile attractive force as 
independent mechanisms. Most recently two contact models were proposed by Thakur 
et al. ED and by Pasha et al. 1421 . which work in the same spirit as Luding’s model, 
but treat loading and un/re-loading behaviors differently. The former excludes the non¬ 
linear elastic stiffness in the plastic regime, and both deal with a more brittle, abrupt 
reduction of the adhesive contact force. The authors further used their models to study 
the scaling and effect of DEM parameters in an uniaxial compression test j43l . and 
compared part of their results with other models l42l . 

When two particles interact, their behavior is intermediate between the extremes 
of perfectly elastic and fully inelastic, possibly even fragmenting, where the latter is 
not considered in this study. Considering a dynamic collision is our choice here, but 
without loss of generality, most of our results can also be applied to a slow, quasi-static 
loading-unloading cycle that activates the plastic loss of energy, by replacing kinetic 

lr The details on the geometry are explained in Refs. [HID. 
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with potential energies. Rozenblat et al. 04) have recently proposed an empirical 
relation between impact velocity and static compression force. 

The amount of energy dissipated during a collision can be best quantified by the 
coefficient of restitution, which is the ratio of magnitude of post-collision and pre¬ 
collision normal relative velocities of the particles. It quantifies the amount of energy 
that is not dissipated during the collision. For the case of plastic and viscoelastic col¬ 
lisions, it was suggested that dissipation depends on impact velocity EM]; this can 
be realized by viscoelastic forces Il46ll48] - [50ll and follows from plastic deformations 
too |5H . 

Early experimental studies 15211531 on adhesive polystyrene latex spheres of mi¬ 
crometer size showed sticking of particles for velocities below a threshold and an in¬ 
creasing coefficient of restitution for velocities increasing above the threshold. Wall 
et al. l54l further confirmed these findings for highly mono-disperse ammonium par¬ 
ticles. Thornton et al. m and Brilliantov et al. EH proposed an adhesive visco- 
elasto-plastic contact model in agreement with these experiments. Work by Sorace 
et al. f56l also confirms the sticking at low velocities for particle sizes of the order 
of a few mm. Li et al. f57l proposed a dynamical model based on JKR for the im¬ 
pact of micro-sized spheres with a flat surface, whereas realisitc particle contacts are 
usually not flat |f58| . Recently, Saitoh et al. |f59l even reported negative coefficients 
of restitution in nanocluster simulations, which is an artefact of the wrong definition 
of the coefficient of restitution; one has to relate the relative velocities to the normal 
directions before and after collision and not just in the frame before collision, which 
is especially a serious effect for softer particles M- Jasevicius et al. 1161116211 have re¬ 
cently studied the rebound behavior of ultrafine silica particles using the contact model 
by Tomas ll22l424l [6?l . They found that energy absorption due to attractive forces 
is the main source of energy dissipation at lower impact velocities or compression, 
while plastic deformation-induced dissipation becomes more important with increas¬ 
ing impact velocity. They found some discrepancies between numerical and experi¬ 
mental observations and concluded that these might be due to the lack of knowledge 
of particle- and contact-parameters, including surface roughness, adsorption layers on 
particle surfaces, and microscopic material property distributions (inhomogeneities), 
which in essence are features of the meso-particles that we aim to study. 

In a more recent study, Shinbrot et al. Il64l studied charged primary particles with 
interesting single particle dynamics in the electromagnetic field. They found ensembles 
of attractive (charged) particles can forming collective contacts or even fingers, extend¬ 
ing the concepts of “contact” well beyond the idealized picture of perfect spheres, as 
shown also in the appendix of the present study. 

Finally, Rathbone et al. f65l presented a new force-displacement law for elasto- 
plastic materials and compare it to their FEM results that resolve the deformations in 
the particle contact zone. This was complemented by an experimental study comparing 
various models and their influence on the bulk flow behavior JT] . 

1.2 Model classification 

Since our main focus is on dry particles, here we do not review the diverse works 
that involve liquid ll66l or strong solid bridges l67l . Even though oblique collisions 
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between two particles are of practical relevance and have been studied in detail by 
Thornton et al. (68ll69) . here we focus on central normal collisions without loss of 
generality. Finally, we also disregard many minute details of non-contact forces, as, 
e.g. due to van der Waals forces, for the sake of brevity, but will propose a very simple 
mesoscale non-contact force model in section 1231 

Based on our review of adhesive, elasto-visco-plastic contact models, here we pro¬ 
pose a possible classification, by dividing them into three groups (based on their com¬ 
plexity and aim): 

(1) Academic contact models, 

(2) Mesoscopic contact models, and 

(3) Realistic, fully detailed contact models. 

Here we focus on adhesive elastic, and elasto-plastic contact models mainly, while the 
effect of various forces on adhesion of fine particles is reviewed in Ref. lITOl . and some 
of the more complex models are reviewed and compared in Ref. [69). 

1. Academic contact models allow for easy analytical solution, as for example 
the linear spring-dashpot model (50) , or piece-wise linear models with constant 
unloading stiffness (see e.g. Walton and Braun (TP ). which feature a constant 
coefficient of restitution (independent of impact velocity). Also the Hertzian 
visco-elastic models belong to this class, even though they provide a velocity 
dependent coefficient of restitution, for a summary see Ref. [50l and references 
therein, while for a recent comparison see Ref. GS- However, no academic 
model can fully describe realistic, practically relevant contacts. Either the ma¬ 
terial or the geometry/mechanics is too idealized; in application, there is hardly 
any contact that is perfectly linear or Hertzian visco-elastic. Academic models 
thus miss most details of real contacts, but can be treated analytically. 

2. Mesoscopic contact models (or, with other words, contact models for meso- 
particles) are a compromise, (i) still rather easy to implement, (ii) aimed for 
fast ensemble/bulk-simulations with many particles and various materials, and 
(iii) contain most relevant mechanisms, but not all the minute details of every 
primary particle and every single contact. They are often piece-wise linear, e.g. 
with a variable unloading stiffness or with an extended adhesive force, leading 
to a variable coefficient of restitution, etc., see Refs. eibsisiieiiiei). 

3. Realistic, full-detail contact models have (i) the most realistic, but often rather 
complicated formulation, (ii) can reproduce with similar precision the pair in¬ 
teraction and the bulk behavior, but (iii) are valid only for the limited class of 
materials they are particularly designed for, since they do include all the minute 
details of these interactions. A few examples include: 

(a) visco-elastic models: Walton (74), Brilliantov (55][75), Haiat 1(761 : 

(b) adhesive elastic models: JKR Ifl6l . Dahneke (771 . DMT IfTT) , Thornton 
and Yin (20); 

(c) adhesive elasto-plastic models: Molerus Q8), Thornton and Ning l2P . 
Tomas l22l424l[63l , Pasha et al. [42l . 
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While the realistic models are designed for a special particulate material in mind, our 
main goal is to define and apply mesoscopic contact models to simulate the bulk be¬ 
havior of a variety of assemblies of many particles (for which no valid realistic model 
is available), we focus on the second class: mesoscopic contact models. 

1.3 Focus and Overview of this study 

In particular, we study the dependence of the coefficient of restitution for two meso- 
particles on impact velocity and contact/material parameters, for a wide range of im¬ 
pact velocities, using the complete version of the contact model by Luding 11131 , with 
a specific piece-wise linear non-contact force term. We observe sticking of parti¬ 
cles at low velocity, which is consistent with previous theoretical and experimental 
works EDEUHE Pasha et al. l42l recently also reproduced the low velocity stick¬ 
ing using an extension of the similar, but simpler model l78l . Above a certain small 
velocity, dissipation is not strong enough to dissipate all relative kinetic energy and 
the coefficient of restitution begins to increase. We want to understand the full regime 
of relative velocities, and thus focus also on the less explored intermediate and high 
velocity regimes, as easily accessible in numerical simulations. In the intermediate 
regime, we observe a decrease in the coefficient of restitution, as observed previously 
for idealized/homogeneous particles H211I55I . however the functional behavior is differ¬ 
ent compared to the predictions by Thornton ED . In Anncndix l2.2.41 we show that this 
property can be tuned by simple modifications to our model. Tanaka et al. [79] have 
recently reported similar results, when simulating the collision of more realistic dust 
aggregates, consisting of many thousands of nanoparticles that interact via the JKR 
model. With further increase in impact velocity, we find a second sticking regime due 
to the non-linearly increasing adhesive and plastic dissipation. For even higher veloci¬ 
ties, the second, intermediate sticking regime is terminated by a second rebound regime 
due to the elastic core that can be specified in the model. Finally, since the physical 
systems under consideration also are viscous in nature, we conclude with some simu¬ 
lations with added viscous damping, which is always added on top of the other model 
ingredients, but sometimes neglected in order to allow for analytical solutions. 

An exemplary application of our model that shows the unexpected high velocity 
sticking and rebound regime (which might not be observed in the case of homoge¬ 
neous granular materials) is, the coating process in cold sprays. In these studies, the 
researchers are interested in analyzing the deposition efficiency of the powder on a 
substrate as a function of the impact velocity. Bonding/coating happens when the im¬ 
pact velocity of the particles exceeds a “critical velocity”, with values of the order of 
10 2 m/s E®. Interestingly, when the velocity is further increased the particles do 
not bond (stick) to the substrate anymore, and a decrease in the deposition efficiency 
(inverse of the coefficient of restitution) is observed 0. Schmidt et al. El have used 
numerical simulations to explore the effect of various material properties on the crit¬ 
ical velocity, while Zhou et al. 0 studied the effect of impact velocity and material 
properties on the coating process, showing that properties of both particle and substrate 
influence the rebound. Using our model, one could explore the dependence of the de¬ 
position efficiency on the impact velocity, leading to the synergy between different 
communities. 
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The paper is arranged as follows: In section [3 we introduce the DEM simulation 
method and the basic contact models for the normal direction; one type of meso-models 
is further elaborated on in the following section [T] where the coefficient of restitution is 
computed analytically, and dimensionless contact parameters are proposed in section^] 
The limit of negligible non-contact forces is considered in section[3 where various spe¬ 
cial cases are discussed, the contact model parameters are studied, and also asymptotic 
solutions and limit values are given, before the study is concluded in section [6] 

2 Discrete Element Method 

The elementary units of particulate systems as granular materials or powders are grains 
that deform under applied stress. Since the realistic and detailed modeling of real par¬ 
ticles in contact is too complicated, it is necessary to relate the interaction force to 
the overlap 8 between two particles in contact. Note that the evaluation of the inter¬ 
particle forces based on the overlap may not be sufficient to account for the inhomo¬ 
geneous stress distribution inside the particles, for internal re-arrangements 126ft . and 
for possible multi-contact effects l45l . However, this price has to be paid in order to 
simulate large samples of particles with a minimal complexity and still taking various 
physical contact properties such as non-linear contact elasticity, plastic deformation or 
load-dependent adhesion into account. 

2.1 Equations of Motion 

If all forces acting on a spherical particle p, either from other particles, from bound¬ 
aries or externally, are known - let their vector sum be f p - then the problem is reduced 
to the integration of Newton’s equations of motion for the translational degrees of free¬ 
dom (the rotational degrees are not considered here since we focus only on normal 
forces) for each particle: m p ^r p = f p + m p g where, m p is the mass of particle p, r p 
its position, f p = YLcfp ' s the total force due to all contacts c, and g is the acceler¬ 
ation due to volume forces like gravity. With tools as nicely described in textbooks 
as lf80U82l . the integration over many time-steps is a straightforward exercise. The 
typically short-ranged interactions in granular media allow for further optimization by 
using linked-cell (LC) or alternative methods in order to make the neighborhood search 
more efficient H831184I . However, such optimization issues are not of concern in this 
study, since only normal pair collisions are considered. 

2.2 Normal Contact Force Laws 

Two spherical particles i and j, with radii a, and aj, r, and rj being the position vectors 
respectively, interact if their overlap, 

8 = (ai + aj)-(ji — rj)-n , (1) 

is either positive, 8 > 0, for mechanical contact, or smaller than a cut-off, 0 > 8 > 8a, 
for non-contact interactions, with the unit vector n = n !; - = (?,■ — rj )/|? ( - — r,| pointing 


r = f^ 



(a) 



Figure 1: Schematic plots of contact forces for (a) the linear normal model for a per¬ 
fectly elastic collision, and (b) the force-overlap relation for an elasto-plastic adhesive 
collision 


from j to i. The force on particle i, from particle j, at contact c, can be decomposed 
into a normal and a tangential part as f c '■= f = f n n + ft, where n ■ t = 0 , n and t 
being normal and tangential parts respectively. In this paper, we focus on frictionless 
particles, i.e., only normal forces will be considered, for tangential forces and torques, 
see e.g. Ref. m and references therein. 

In the following, we discuss various normal contact force models, as shown schemat¬ 
ically in Fig. 1. We start with the linear contact model (Fig. 1 (a)) for non-adhesive par¬ 
ticles, before we introduce a more complex contact model that is able to describe the 
realistic interaction between adhesive, inhomogeneous 0 , slightly non-spherical parti¬ 
cles (Fig. 1(b)). 

2.2.1 Linear Normal Contact Model 

Modelling a force that leads to an inelastic collision requires at least two ingredients: 
repulsion and some sort of dissipation. The simplest (but academic) normal force law 
with the desired properties is the damped harmonic oscillator 

f n =k8 + yov" , ( 2 ) 

with spring stiffness k, viscous damping yo, and normal relative velocity v” = — v,j • n = 
— (vj — Vj) -n = 8. This model (also called linear spring dashpot (LSD) model) has the 
advantage that its analytical solution (with initial conditions 8 (0) =0 and 8 (0 ) = Vq) 
allows easy calculations of important quantities ll50l . For the non-viscous case, the 
linear normal contact model is given schematically in Fig. 1. 

2 Examples of inhomogeneous particles include core-shell materials, clusters of fine primary particles or 
randomly micro-porous particles. 
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The typical response time (contact duration) and the eigenfrequency of the contact 
are related as 


t c = ^ and co = ^ ( k/m r ) - t] ( j 


(3) 


with the rescaled damping coefficient rjo = yo/(2 m r ), and the reduced mass m r = 
nijirij/(nij + nij), where the rjo is defined such that it has the same units as co, i.e., 
frequency. From the solution of the equation of a half-period of the oscillation, one 
also obtains the coefficient of restitution 


e « S ° = V /M = exp (- 7 : 770 /®) = exp(-i)ofc) , (4) 


which quantifies the ratio of normal relative velocities after (v/) and before (v,) the 
collision. Note that in this model e„ is independent of v,-. For a more detailed review 
on this and other, more realistic, non-linear contact models, see mm and references 
therein. 

The contact duration in Eq. d3} is also of practical and technical importance, since 
the integration of the equations of motion is stable only if the integration time-step At 
is much smaller than t c . Note that t c depends on the magnitude of dissipation: In the 
extreme case of an over-damped spring (high dissipation), t c can become very large 
(which renders the contact behavior artificial l48l ). Therefore, the use of neither too 
weak nor too strong viscous dissipation is recommended, so that some artificial effects 
are not observed by the use of viscous damping. 


2.2.2 Adhesive Elasto-PIastic Contacts 

For completeness, we re-introduce the piece-wise linear hysteretic model lfl5l as an 
alternative to non-linear spring-dashpot models or more complex hysteretic models 
ll2T1j24l[85l[86l . It reflects permanent plastic deformation 0, which takes place at the 
contact, and the non-linear increase of both elastic stiffness and attractive (adhesive) 
forces with the maximal compression force. 

In Fig. 2, the normal force at contact is plotted against the overlap 8 between two 
particles. The force law can be written as 

( k\8 if k 2 (8 — 8o)>k\8 

/ hys = < k 2 (8-8 0 ) if kiS >k 2 (S-S 0 ) > -k c S (5) 

[ -k c 8 if -k c 8 > k 2 (8- 8 q) 

with k\ <k 2 < k p , respectively the initial loading stiffness, the un-/re-loading stiffness 
and the elastic limit stiffness. The latter defines the limit force branch k p {8 — 8q), 
as will be motivated next in more detail, and k 2 interpolates between k\ and k p , see 
Eq. ©. For k, = 0, the above contact model reduces to that proposed by Walton and 
Braun l ITTl , with coefficient of restitution 

eT = Vhjh • ( 6 ) 

3 After a contact is opened, the pair forgets its previous contact, since we assume that the contact points 
at a future re-contact of the same two particles are not the same anymore. 
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Figure 2: Schematic graph of the piece-wise linear, hysteretic, and adhesive force- 
displacement model in normal direction from Ref. m. 
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During the initial loading the force increases linearly with overlap 8 along k \, until 
the maximum overlap <5 max = v< \Jm r /k\ (for binary collisions) is reached, which is a 
history parameter for each contact. During unloading the force decreases along k 2 , see 
Eq. ([9]), from its maximum value k\ <5 max at 5 max down to zero at overlap 

8 0 = (l-k l /k 2 )8 maLK , (7) 


where So resembles the permanent plastic contact deformation. Re-loading at any 
instant leads to an increase of the force along the (elastic) branch with slope k 2 , until the 
maximum overlap <5 max (which was stored in memory) is reached; for still increasing 
overlap 8, the force again increases with slope k\ and the history parameter 5 max has to 
be updated. 

Unloading below 8q leads to a negative, attractive (adhesive) force, which follows 
the line with slope k 2 , until the extreme adhesive force —k c 8 m ; n is reached. The corre¬ 
sponding overlap is 


_ (k 2 -k 1 ) 
{k 2 + k c ) 


( 8 ) 


Further unloading follows the irreversible tensile limit branch, with slope — k c , with the 
attractive force / hys = —k c S. 

The lines with slope k\ and —k c define the range of possible positive and negative 
forces. Between these two extremes, unloading and/or re-loading follow the line with 
slope k 2 . A non-linear un-/re-loading behavior would be more realistic, however, due 
to a lack of detailed experimental informations, the piece-wise linear model is used as 
a compromise, besides that it is easier to implement. The elastic k 2 branch becomes 
non-linear and ellipsoidal if a moderate normal viscous damping force is active at the 
contact, as in the LSD model. 

In order to account for realistic load-dependent contact behavior, the k 2 value is 
chosen to depend on the maximum overlap <5 max , i.e. contacts are stijfer and more 
strongly plastically deformed for larger previous deformations so that the dissipation 
depends on the previous deformation history. The dependence of k 2 on overlap <5 max is 
chosen empirically as linear interpolation: 


k 2 (&max) — 


kp if <5 max /<S£ ax > 1 

k\ + {kp k 1 ) 8 max l ^m ax 

if <5 ma x/$£ ax < 1 


(9) 


where k p is the maximal (elastic) stiffness, and 


sp kp 2 a\a 2 

<?m ax = , , 0/ , 

k v -k 1 a\+a 2 


( 10 ) 


is the plastic flow limit overlap, with (j)f the dimensionless plasticity depth, a\ and a 2 
being the radii of the two particles. This can be further simplified to 


t>o = 0/ai2, 


( 11 ) 


where 8 ( \' represents the plastic contact deformation at the limit overlap, and a\ 2 = 
is the reduced radius. In the range <5 max < t>m ax , the stiffness k 2 can also be 
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written as: 


k 2 =k ] + ik ‘ , J' ) n\ ( 12 ) 

k i o max 

where / max = k\ 5 max is the same as Eq. (4) in fTTI with prefactor S = /* - 1 . 

M ^max 

From energy balance considerations, one can define the “plastic” limit velocity 

v p = k\/m r <5,£ ax , (13) 

below which the contact behavior is elasto-plastic, and above which the perfectly elas¬ 
tic limit-branch is reached. Impact velocities larger than v p have consequences, as 
discussed next (see Sec. 12.2.41 ). 

In summary, the adhesive, elasto-plastic, hysteretic normal contact model is defined 
by the four parameters k\, k p , k c and <j)f that, respectively, account for the initial plastic 
loading stiffness, the maximal, plastic limit (elastic) stiffness, the adhesion strength, 
and the plastic overlap-range of the model. This involves an empirical choice for the 
load-dependent, intermediate elastic branch stiffness kz, which renders the model non¬ 
linear in its behavior (i.e. higher confinging stress leads to stiffer contacts like in the 
Hertz model), even though the present model is piece-wise linear. 

2.2.3 Motivation of the original contact model 

To study a collision between two ideal, homogeneous spheres, one should refer to 
realistic, full-detail contact models with a solid experimental and theoretical foundation 
queued. These contact models feature a small elastic regime and the particles 
increasingly deform plastically with increasing, not too large deformation (overlap). 
During unloading, their contacts end at finite overlap due to flattening. An alternative 
model was recently proposed, see Ref. l42l . that follows the philosophy of plastically 
flattened contacts with instantaneous detachment at positive overlaps. 

However, one has to also consider the possibility of rougher contacts [58 j, and 
possible non-contact forces that are usually neglected for very large particles, but can 
become dominant and hysteretic as well as long-ranged for rather small spheres ll22l 

M- 

The mesoscopic contact model used here, as originally developed for sintering 
j25l , and later defined in a temperature-independent form 031 , follows a different ap¬ 
proach in two respects: (i) it introduces a limit to the plastic deformation of the par¬ 
ticles/material for various reasons as summarized below and also in subsection 12.2.41 
and (ii) the contacts are not idealized as perfectly flat, and thus do not have to lose 
mechanical contact immediately at un-loading, as will be detailed in subsection l2.2.5l 
Note that a limit to the slope k p resembles a simplification of different contact be¬ 
havior at large deformations'. 

(i) for low compression, due to the wide probability distribution of forces in bulk gran¬ 
ular matter, only few contacts should reach the limit, which would not affect much the 
collective, bulk behavior; 

(ii) for strong compression, in many particle systems, i.e., for large deformations, the 
particles cannot be assumed to be spherical anymore, and they deform plastically or 
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could even break; 

(iii) from the macroscopic point of view, too large deformations would lead to volume 
fractions larger than unity, which for most materials (except highly micro-porous, frac¬ 
tal ones) would be unaccountable; 

(iv) at small deformation, contacts are due to surface roughness realized by multi¬ 
ple surface asperities and at large deformation, the single pair point-contact argument 
breaks down and multiple contacts of a single particle can not be assumed to be inde¬ 
pendent anymore; 

(v) finally, (larger) meso-particles have a lower stiffness than (smaller) primary parti¬ 
cles ED, which is also numerically relevant, since the time step has to be chosen such 
that it is well below the minimal contact duration of all the contacts. If ki is not limited 
the time-step could become prohibitively small, only because of a few extreme (large 
compression) contact situations. 

The following two subsections discuss the two major differences of the present piece- 
wise linear (yet non-linear) model as compared to other existing models: (i) the elastic 
limit branch, and (ii) the elastic re-loading or non-contact-loss, as well as their reasons, 
relevance and possible changes/tuning - in case needed. 

2.2.4 Shortcomings, physical relevance and possible tuning 

In the context of collisions between perfect homogeneous elasto-plastic spheres, a 
purely elastic threshold/limit and enduring elastic behavior after a sharply defined 
contact-loss are indeed questionable, as the plastic deformation of the single particle 
cannot become reversible/elastic. Nevertheless, there are many materials that support 
the idea of a more elastic behavior at large compression (due to either very high im¬ 
pact velocity or multiple strong contact forces), as discussed further in the paragraphs 
below. 

Mesoscopic contact model applied to real materials: First we want to recall that 
the present model is mainly aimed to reproduce the behavior of multi-particle systems 
of realistic fine and ultra-fine powders, which are typically non-spherical and often 
mesoscopic in size with internal micro-structure and micro-porosity on the scale of 
typical contact deformation. For example, think of clusters/agglomerates of primary 
nano-particles that form fine micron-sized secondary powder particles, or other fluffy 
materials l26l . The primary particles are possibly better described by other contact 
models, but in order to simulate a reasonable number of secondary (meso) particles 
one cannot rely on this bottom-up approach and hence a mesoscopic contact model 
needs to be used. During the bulk compression of such a system, the material deforms 
plastically and both the bulk and particles’ internal porosity reduces j26l . Plastic defor¬ 
mation diminishes if the material becomes so dense, with minimal porosity, such that 
the elastic/stiff primary particles dominate. Beyond this point the system deforms more 
elastically, i.e. the stiffness becomes high and the (irrecoverable) plastic deformations 
are much smaller than at weaker compression. 

In their compression experiments of granular beds with micrometer sized granules 
of micro-crystalline cellulose, Persson et al. l87l found that a contact model where a 
limit on plastic deformation is introduced can very well describe the bulk behavior. 
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Experimentally they observe a strong elasto-plastic bulk-behavior for the assembly at 
low compression strain/stress. In this phase the height of the bed decreases, irreversibly 
with the applied load. It becomes strongly non-linear beyond a certain strain/stress, 
which is accompanied by a dramatic increase of the stiffness of the aggregate. They 
associate this change in the behavior to the loss of porosity and the subsequent more 
elastic bulk response to the particles that are now closely in touch with each other. In 
this new, re-structured, very compacted configurations, further void reduction is not 
allowed anymore and thus the behavior gets more elastic. While the elastic limit in the 
contact model does not affect the description of the bulk behavior in the first part, the 
threshold is found to play a key role in order to reproduce the material stiffening (see 
Fig. 8 in Ref. 11871 ). 

Note that in an assembly of particles, not all the contacts will reach the limit branch 
and deform elastically simultaneously. That is, even if few contacts are in the elastic 
limit, the system will always retain some plasticity, hence the assembly will never be 
fully elastic. 

Application to pair interactions: Interestingly, the contact model in Sec. 12.2.21 is 
suitable to describe the collision between pairs of particles, when special classes of 
materials are considered, such that the behavior at high velocity and thus large defor¬ 
mation drastically changes. 

(i) Core-shell materials. The model is perfectly suited for plastic core-shell materi¬ 
als, such as asphalt or ice particles, having a “soft” plastic outer shell and a rather stiff, 
elastic inner core. For such materials the stiffness increases with the load due to an in¬ 
creasing contact surface. For higher deformations, the inner cores can come in contact, 
which turns out to be almost elastic when compared to the behavior of the external 
shell. The model was successfully applied to model asphalt, where the elastic inner 
core is surrounded by a plastic oil or bitumen layer {88] ■ Alternatively, the plastic shell 
can be seen as the range of overlaps, where the surface roughness and inhomogeneities 
lead to a different contact mechanics as for the more homogeneous inner core. 

(ii) Cold spray. An other interesting system that can be effectively reproduced by 
introducing an elastic limit in the contact model is cold spray. Researchers have ex¬ 
perimentally and numerically shown that spray-particles rebound from the substrate at 
low velocities, while they stick at intermediate impact energy Il2l44l [89l . Wu et al. a 
experimentally found that rebound re-appears with a further increase in velocity (Fig. 
3 in Ref. 0). Schmidt et al. 0 relate the decrease of the deposition efficiency (in¬ 
verse of coefficient of restitution) to a transition from a plastic impact to hydrodynamic 
penetration (Fig. 16 in Ref. 0). Recently, Moridi et al. Il89l numerically studied the 
sticking and rebound processes, by using the adhesive elasto-plastic contact model of 
Fuding on, and their prediction of the velocity dependent behavior is in good agree¬ 
ment with experiments. 

(iii) Sintering. As an additional example, we want to recall that the present meso¬ 
scopic contact model has already been applied to the case of sintering, see Ref. ]25ll90l . 
For large deformations, large stresses, or high temperatures, the material goes to a 
fluid-like state rather than being solid. Hence, the elasticity of the system (nearly in¬ 
compressible melt) determines its limit stiffness, while tj)f determines the maximal 


15 






volume fraction that can be reached. 

All the realistic situations described above clearly hint at a modification in the con¬ 
tact phenomenology that can not be described solely by an elasto-plastic model beyond 
some threshold in the overlap/force. The limit stiffness k p and the plastic layer depth 
(j)f in our model allow the transition of the material to a new state. Dissipation on the 
limit branch - which otherwise would be perfectly elastic - can be taken care of, by 
adding a viscous damping force (as the simplest option). Due to viscous damping, the 
unloading and re-loading will follow different paths, so that the collision will never be 
perfectly elastic, which is in agreement with the description in Jasevicius et al. mm 
and will be shown later in AppendixlBl 

Finally, note that an elastic limit branch is surely not the ultimate solution, but 
a simple first model attempt - possibly requiring material- and problem-adapted im¬ 
provements in the future. 

Tuning of the contact model: The change in behavior at large contact deformations 
is thus a feature of the contact model which allows us to describe many special types 
of materials. Nevertheless, if desired (without changing the model), the parameters can 
be tuned in order to reproduce the behavior of materials where the plasticity increases 
with deformation without limits, i.e., the elastic core feature can be removed. The 
limit-branch where plastic deformation ends is defined by the dimensionless parame¬ 
ters plasticity depth, (j)f , and maximal (elastic) stiffness, k p . Owing to the flexibility of 
the model, it can be tuned such that the limit overlap is set to a much higher value that 
is never reached by the contacts. When the new value of (j)f is chosen, a new k p can 
be calculated to describe the behavior at higher overlap (as detailed in Appendix^. In 
this way the model with the extended (j)f exhibits elasto-plastic behavior for a higher 
velocity/compression-force range, while keeping the physics of the system for smaller 
overlap intact. 

2.2.5 Irreversibility of the tensile branch 

Finally we discuss a feature of the contact model in G3 , that postulates the irreversibil¬ 
ity, i.e. partial elasticity, of the tensile k c branch, as discussed in Sec. 12.2.21 While this 
is unphysical in some situations, e.g. for homogeneous plastic spheres, we once again 
emphasize that we are interested in non-homogeneous, non-spherical meso-particles, 
as e.g. clusters/agglomerates of primary particles in contact with internal structures of 
the order of typical contact deformation. The perfectly flat surface detachment due to 
plasticity happens only in the case of ideal, elasto-plastic adhesive, perfectly spherical 
particles (which experience a large enough tensile force). In almost all other cases, the 
shape of the detaching surfaces and the hence the subsequent unloading behavior de¬ 
pends on the relative strengths of plastic dissipation, attractive forces, and various other 
contact mechanisms. In the case of meso-particles such as the core-shell materials ||88l , 
assemblies of micro-porous fine powders 112611871 . or atomic nanoparticles ED, other 
details such as rotations can be important. We first briefly discuss the case of ideal 
elasto-plastic adhesive particles and later describe the behavior of many particle sys¬ 
tems, which is the main focus of this work. 
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Ideal homogeneous millimeter sized particles detach with a permanently flattened 
surface created during deformation are well described using contact models presented 
in HUM- This flattened surface is of the order of micrometers and the plastic dis¬ 
sipation during mechanical contact is dominant over the van der Waals force. During 
unloading, when the particles detach, the force suddenly drops to zero from the ten¬ 
sile branch. When there is no contact, further un- and re-loading involves no force. 
Even when the contact is re-established, the contact is still assumed to be elastic, i.e., it 
follows the previous contact-unloading path. This leads to very little or practically no 
plastic deformation at the re-established contact, until the (previously reached) maxi¬ 
mum overlap is reached again and the plasticity kicks in. 

On the other hand for ultra-fine ideal spherical particles of the order of macro¬ 
meters 122116311911 . the van der Waals force is much stronger and unloading adhesion 
is due to purely non-contact forces. Therefore, the non-contact forces do not vanish 
and even extend beyond the mechanical first contact distance. The contact model of 
Tomas 112211631 is reversible for non-contact and features a strong plastic deformation 
for the re-established contact - in contrast to the previous case of large particles. 

The contact model by Luding 113 follows similar considerations as others, ex¬ 
cept for the fact that the mechanical contact does not detach (for details see the next 
section). The irreversible, elastic re-loading before complete detachment can be seen 
as a compromise between small and large particle mechanics, i.e. between weak and 
strong attractive forces. It also could be interpreted as a premature re-establishment 
of mechanical contact, e.g. due to a rotation of the deformed, non-spherical particles. 
Detachment and remaining non-contact is only then valid if the particles do not rotate 
relative to each other; in case of rotations, both sliding and rolling degrees of freedom 
can lead to a mechanical contact much earlier than in the ideal case of a perfect normal 
collision of ideal particles. In the spirit of a mesoscopic model, the irreversible contact 
model is due to the ensemble of possible contacts, where some behave as imagined 
in the ideal case, whereas some behave strongly different, e.g. due to relative rotation. 
However, there are several other reasons to consider an irreversible unloading branch, 
as summarized in the following. 

In the case of asphalt (core-shell material with a stone core and bitumen-shell), 
depending on the composition of the bitumen (outer shell), which can contain a con¬ 
siderable amount of fine solid, when the outer shells collide the collision is plastic. In 
contrast, the collision between the inner cores is rather elastic (even though the inner 
cores collide when the contact deformation is very large). Hence, such a material will 
behave softly for loading, but will be rather stiff for re-loading (elastic k .2 branch), since 
the cores can then be in contact. A more detailed study of this class of materials goes 
beyond the scope of this study and the interested reader is referred to Ref. l88l 

For atomistic nano-particles and for porous materials, one thing in common is the 
fact that the scale of a typical deformation can be much larger than the inhomogeneities 
of the particles and that the adhesion between primary particles is strong enough to 
keep them agglomerated during their re-arrangements (see Fig. 5 in Ref. fZ2) and 
the phenomenology in Ref. l26l . as well as recent results for different deformation 
modes lf27l ). Thus the deformation of the bulk material will be plastic (irreversible), 
even if the primary particles would be perfectly elastic. 

For agglomerates or other mesoscopic particles, we can not assume permanent ideal 
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flattening and complete, instantaneous loss of mechanical contact during unloading 
j26l . In average, many contacts between particles might be lost, but - due to their 
strong attraction - many others will still remain in contact. Strong clusters of primary 
particles will remain intact and can form threads, a bridge or clumps during unloading 
- which either keeps the two surfaces in contact beyond the (idealized) detachment 
point ll26l or can even lead to an additional elastic repulsion due to a clump of particles 
between the surfaces (see Fig. 3 in Ref. fl5l and Appendix[F]). 

During re-loading, the (elastic) connecting elements influence the bulk response. At 
the same time, the re-arrangements of the primary particles (and clusters) can happen 
both inside and on the surface, which leads to reshaping, very likely leaving a non-flat 
contact surface lfTl l26l[58l . As often mentioned for granular systems, the interaction 
of several elastic particles does not imply bulk elasticity of the granular assembly, due 
to (irreversible) re-arrangements in the bulk material - especially under reversal of 
direction [35j. Thus, in the present model an irreversible tensile branch is assumed, 
without distinction between the behavior before and after the first contact-loss-point 
other than the intrinsic non-linearity in the model: The elastic stiffness for re-loading 
k 2 decreases the closer it comes to 8 = 0; in the present version of the contact model, 
k 2 for unloading from the k\ branch and for re-loading from the k c branch are exactly 
matched (for the sake of simplicity). 

It is also important to mention that large deformation, and hence large forces are 
rare, thanks to the exponential distribution of the deformation and thus forces, as shown 
by our studies using this contact model Il25ii29ll92l l. Hence, such large deformations 
are rare and do not strongly affect the bulk behavior, as long as compression is not too 
strong. 

As a final remark, for almost all models on the market - due to convenience and 
numerical simplicity, in case of complete detachment 8 < 0 - the contact is set to its 
initial state, since it is very unlikely that the two particles will touch again at exactly 
the same contact point as before. On the other hand in the present model a long- 
range interaction is introduced, in the same spirit as 02311631 . which could be used to 
extend the contact memory to much larger separation distances. Re-loading from a 
non-contact situation (5 < 0) is, however, assumed to be starting from a “new” contact, 
since contact model and non-contact forces are considered as distinct mechanisms, for 
the sake of simplicity. Non-contact forces will be detailed in the next subsection. 

2.3 Non-contact normal force 

It has been shown in many studies that long-range interactions are present when dry 
adhesive particles collide, i.e. non-contact forces are present for negative overlap 8 
IU51l21l[63ll93ll94 l. In the previous section, we have studied the force laws for contact 
overlap 8 > 0. In this section we introduce a description for non-contact, long range, 
adhesive forces, focusing on the two non-contact models schematically shown in Fig. 
3 - both piece-wise linear in the spirit of the mesoscopic model - namely the reversible 
model and the jump-in (irreversible) non-contact models (where the latter could be seen 
as an idealized, mesoscopic representation of a liquid bridge, just for completeness). 
Later, in the next section, we will combine non-contact and contact forces. 
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Figure 3: Schematic plots of (a) the non-contact adhesive force-overlap relation and 
(b) the non-contact jump-in force-overlap relation. 

2.3.1 Reversible Adhesive force 

In Fig. 3(a) we consider the reversible attractive case, where a (linear) van der Waals 
type long-range adhesive force is assumed. The force law can be written as 



if 5>0 
if 0 > 8 > d a 
if 8 a > 8 


(14) 


with the range of interaction 8 a = —f a /K- < 0 , where k a c > 0 is the adhesive “stiff¬ 
ness” of the material 0 and f a > 0 is the (constant) adhesive force magnitude, active 
also for overlap 8 > 0, in addition to the contact force. When 5 = 0 the force is —f a . 
The adhesive force / adh is active when particles are closer than 8„, when it starts to in¬ 
crease/decrease linearly along — k%, for approach/separation, respectively. In the results 
and theory part of the paper, for the sake of simplicity and without loss of generality, 
the adhesive stiffness can be either chosen as infinite, which corresponds to zero range 
non-contact force ( 8 a = 0 ), or as coincident with the contact adhesive stiffness, e.g. in 
Sec. 12.2.21 that is k° = k c . 

2.3.2 Jump-in (Irreversible) Adhesive force 

In Fig. 3(b) we report the behavior of the non-contact force versus overlap when the 
approach between particles is described by a discontinuous (irreversible) attractive law. 
The jump-in force can be simply written as 



if 5 <0 
fa if 5 > 0 


(15) 


As suggested in previous studies Ifl6ll21ll55l . there is no attractive force before the 
particles come into contact; the adhesive force becomes active and suddenly drops to 

4 Since the i c -branch has a negative slope, this parameter does not represent a true stiffness of the material, 
which must have a positive sign. 
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a negative value, —f a , at contact, when 8 = 0. The jump-in force resembles the limit 
case —> °° of Eq. Note that the behavior is defined here only for approach of 

the particles. We assume the model to be irreversible, as in the unloading stage, during 
separation, the particles will not follow this same path (details will be discussed below). 

3 Coefficient of Restitution 

The amount of dissipated energy relative to the incident kinetic energy is quantified by 
1 — e 2 , in terms of the coefficient of restitution e. Considering a pair collision, with 
particles approaching from infinite distance, the coefficient of restitution is defined as 

e = -5— (16a) 

v,-" 

where vf° and v,°° are final and initial velocities, respectively, at infinite separations 
(distance beyond which there is no long range interaction). Assuming superposition 
of the non-contact and contact forces, the restitution coefficient can be further decom¬ 
posed including terms of final and initial velocities, Vf and v<, at overlap 5 = 0, where 
the mechanical contact-force becomes active: 

Vf°° Vf Vf 

e = — - J -— = e 0 e n £i, (16b) 

Vf Vi V( 

and £j and £„ are the pull-in and pull-off coefficients of restitution, that describe the 
non-contact parts of the interaction (8 < 0), for approach and separation of particles, 
respectively. The coefficient of restitution for particles in mechanical contact (5 >0) 
is e„, as analytically computed in subsection l3.3l 

In the following, we will analyze each term in Eq. (1 1 6bb separately, based on energy 
considerations. This provides the coefficient of restitution for a wide, general class 
of meso interaction models with superposed non-contact and contact components, as 
defined in sections I2?2ll2.31 

For the middle term, e n , different contact models with their respective coefficients 
of restitution can be used, e.g. eJ; SD from Eq. (Q}, e^ B from Eq. or e)/ YS as calcu¬ 
lated below in subsec tion !3. 31 Prior to this, we specify e, in subsection 13 T I and then £ a 
in subsection l3.2l for the simplest piece-wise linear non-contact models. |f| 


3.1 Pull-in coefficient of restitution 


In order to describe the pull-in coefficient of restitution £,-, we focus on the two non- 
contact models proposed in Sec. 12.31 as simple interpretations of the adhesive force 
during the approach of the particles. 

When the reversible adhesive contact model is used, energy conservation leads to 
an increase in velocity due to the attractive branch from 8 U (< 0) to contact: 


1 

-m r v, 


a 


-J a 8 a - 


1 


-m r Vi 


(17a) 


5 If other, possibly non-linear non-contact forces such as square-well, van der Waals or Coulomb are 
used, see Refs. (SHIS, the respective coefficient of restitution has to be computed, and also the long-range 
nature has to be accounted for, which goes far beyond the scope of this paper. 
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which yields 


„rev—adh 


Vj 




fl/K 


(17b) 


The pull-in coefficient of restitution is thus larger than unity; it increases with increas¬ 
ing adhesive force magnitude f a and decreases with the adhesive strength of the mate¬ 
rial k“ (which leads to a smaller cutoff distance). 

On the other hand, if the irreversible adhesive jump-in model is implemented, a 
constant value £- ump m = 1 is obtained for first approach of two particles, before con¬ 
tact, as yJ um P- ln — o for 8 < 0 and the velocity remains constant v; = v, 00 . 


3.2 Pull-off coefficient of restitution 


The pull-off coefficient of restitution is defined for particles that lose contact and sepa¬ 
rate. Using the adhesive reversible model, as described in section [2.3. 1 1 energy balance 
leads to a reduction in velocity during separation: 



-Ja8 a 


IrVf 


(18a) 


which yields 





(18b) 


due to the negative overlap 8 a at which the contact ends. Similar to Eq. ( 1 1 7bb . the 
pull-off coefficient of restitution depends on both the adhesive force magnitude f a and 
stiffness k c , given the separation velocity vj at the end of the mechanical contact. 

It is worthwhile to note that the force-overlap picture described above, with e a < 1 
defined as in Eq. (1 1 8bb refers to a system with sufficiently high impact velocity, so that 
the particles can separate with a finite kinetic energy at the end of collision, i.e., v f 2 > 
fci/( m rK) =: ( v /) 2 or > equivalently, > vjj (e«£,‘), where v a j denotes the minimal 
relative velocity at the end of the contact, for which particles can still separate. If 
the kinetic energy reaches zero before the separation, e.g. the particles start re-loading 
along the adhesive branch until the value 8 = 0 is reached and the contact model kicks 
in. 


3.3 Elasto-plastic coefficient of restitution 


The key result of this paper is the analytical study of the coefficient of restitution as 
function of the impact velocity, for the model presented in Fig. 4(b), disregarding vis¬ 
cous forces in order to allow for a closed analytical treatment. The impact velocity v, is 
considered for two cases v, < v p and v, > v p , with the plastic-limit velocity v p (needed 
to reach the elastic branch), defined as: 


v 


p 



(19) 
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(a) (b) 


Figure 4: (a) Reversible and irreversible non-contact forces, where the top blue line (for 
negative overlap) represents the former and the bottom red line (for negative overlap) 
the latter. The black line for positive overlap represents the linear contact force as 
superimposed on the non-contact force, (b) Force-displacement law for elasto-plastic, 
adhesive contacts superimposed on the irreversible non-contact adhesive force. 


where the term(s) with f a represent the energy gained or lost by this (attractive, neg¬ 
ative) constant force, with zero reached at overlap 5j 11 = f a /k\, and 5max defined in 
Eq. dTOl) . The velocity v p needed to reach the limit branch thus decays with increasing 
non-contact attraction force f a . 


3.3.1 Plastic contact with initial relative velocity v, < v p 

When Vi < v p the particles after loading to 5 max , unload with slope and the system 
deforms along the path 0 — y <5 max — > 8{j — > 5 mm —> 0, corresponding to A —> B —> C — > 
D — » E in Fig. 4(b). 

The initial kinetic energy (at 8 = 0 overlap, with adhesive force f a and with initial 
velocity v, < v p ) is completely transformed to potential energy at the maximum overlap 
<5 ma x where energy-balance provides: 

E t := l-m r vj = f a ) ( 5 max ~ = \^x(k\8 max -2fa) , (20a) 

2 2 \ k\ J 2 k\ 2 


so that the physical (positive) solution yields: 


^max — 


fa + J fa + kim r vj I / (1 \\ 2 — 

-£- = 8a + y i^8a +m r vj/ki , (20b) 


with zero force during loading at 8^ 1 = f a /k\. The relative velocity is reversed at 
<5 m ax. and unloading proceeds from point B along the slope k 2 = AofSmax )■ Part of the 
potential energy is dissipated, the rest is converted to kinetic energy at point C, the 
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force-free overlap 8f in the presence of the attractive force f a : 


— ^(^l^max /a)(<Vax 8q) — (^ n r V i + ^ ^ ^ max / Q ) > 

(20c) 

(Y\ (Y\ 

where 5g = [(k 2 - fcptV ax + f a ]/k 2 =: 8 0 + <5„, with Sy = f a /k 2 , and the second 
identity follows fromEq. (I20ab . using the force balance at the point of reversal k 2 (8 miix — 

<?(?) = (Vax-/«■ 0 

Further unloading, below 8f leads to attractive forces. The kinetic energy at 8{[ is 
partly converted to potential energy at point D, with overlap 5 mm , where the minimal 
(maximally attractive) force is reached. Energy balance provides: 

— Wirt'min = ~^ m r v 0 — ^^ 2(^0 — ^min)” = 7^ m r v Q ~ ^j~{kc-8min + fa) 1 (20d) 

where 5 m ; n = k2 ^° + /“ = ^ 2 ^j_] i: 5max , and the second identity follows from inserting 

8 0 = (1 +^eA2)^min +fa/k 2 . 

The total energy is finally converted exclusively to kinetic energy at point E, the 
end of the collision cycle (with overlap 5 = 0): 


T'" rV / = \ ,llrV min - \ k c 5 ^ n - faS m 


( 20 e) 


Using Eqs. (120cl) . (120dl) . and (120el) with the definition of <S mm , and combining terms 
proportional to powers of f a and t) max yields the final kinetic energy after contact: 


?W.- 


k\ k c (k 2 -k x ) 2 


E f := —m r v~f = -- 

f 2 f [k 2 k\k 2 ( k 2 + k c ) \ 2 


0 ^l<Vax fa8 n 


( 21 ) 


with <5 max as defined in Eq. (I20bb . Note that the quadratic terms proportional to f 2 have 
cancelled each other, and that the special cases of non-cohesive ( k c = 0 and/or f a = 0 ) 
are simple to obtain from this analytical form. Finally, dividing the final by the initial 
kinetic energy, Eq. (I20ab . we have expressed the coefficient of restitution 





( 22 ) 


as a function of maximal overlap reached, <5 max . non-contact adhesive force, f a , elastic 
unloading stiffness, k 2 = k 2 (8 nax), and the constants plastic stiffness, k\ , and cohesive 
“stiffness”, k c . 


3.3.2 Plastic-elastic contact with initial relative velocity v,- > v p 

When the initial relative velocity v,- is large enough such that v i > Vp, the estimated 
maximum overlap <5 max as defined in Eq. < 120bl ) is greater than (5£ ax . Let iq be the 

6 From this point, we can derive the coefficient of restitution for the special case of k c = 0 final energy, 
using the final energy Ej}\k c = 0) := (f + 2k\ F,j) /{2ki). 
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velocity at overlap 8 ^ ax . The system deforms along the path 0 —> 8 ^ ax —> 5 max —> 
8 q —> 5 m in —> 0. 

The initial kinetic energy (at 8 = 0 overlap, with adhesive force f a and with initial 
velocity v, > v p ) is not completely converted to potential energy at 8 = 5m ax , where 
energy balance provides: 


\ m A = ^m r vj - ] -m r v 2 p = l -m r v 2 - ^5£ ax (k| 5£ ax - 2f a ) , 


(23a) 


using the definition of v p in Eq. ( IT9l >. 

From this point the loading continues along the elastic limit branch with slope k p 
until all kinetic energy is transferred to potential energy at overlap 5 max > 5,n ax , where 
the relative velocity changes sign, i.e., the contact starts to unload with slope k p . Since 
there is no energy disspated on the k^-branch (in the absence of viscosity), the potential 
energy is completely converted to kinetic energy at the force-free overlap 8q\ on the 
plastic limit branch 


\ m A = (MLx ~fa ) 2 + \™ r v T , 


(23b) 


with the first term taken from Eq. ( I20cb , but replacing <S max with 5, aax and ki by k p . 

Further unloading, still with slope k p , leads to attractive forces. The kinetic energy 
at Sq P is partly converted to potential energy at 8 P mp where energy balance yields: 


Vn r vL a = n r vl l -k p {8? - 8^) 2 = l -m r v 2 - - f a ) 


(23c) 


Some of the remaining potential energy is converted to kinetic energy so that at the 
end of collision cycle (with overlap 5 = 0) one has 

\ m A = ^ k c (5^ in ) 2 - f a 8 p mn , (23d) 

analogously to Eq. (120eb 

When Eq. (123dl) is combined with Eqs. (123bb and (I23cb . and inserting the definitions 

5 min = k, ’k° p +k!“ = %+lf aX ’ and S o P = 0 + k c/k p )8^ ln +fa/k p , one obtains (similar 
to the previous subsection): 


e 1 2) = 


1 


1 


f = ? ,Tl r V f = j m r V i 


k\ | K (k p k\) 2 k\ 2 
kphk p (k p +k c )\ 2 


(23e) 


Dividing the final by the initial kinetic energy, we obtain the coefficient of restitution 


e [ n ] = yjEf/Ei =: y/l -E diss /Ei , (24) 

with constants k\, k p , k c , f a , and 5£ ax . Note that el,~\ interestingly, does not depend on 
f a at all, since the constant energy E &; ss is lost exclusively in the hysteretic loop (not 
affected by f a ). Thus, even though /-’ ( | lxs does not depend on the impact velocity, the 
coefficient of restitution does, because of its definition. 
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As final note, when the elastic limit regime is not used, or modified towards larger 
5max, as defined in appendix [A] the limit velocity, v p , increases, and the energy lost, 
£diss, increases as well (faster than linear), so that the coefficient of restitution just 

( 2 ) 

becomes e„ = 0, due to complete loss of the initial kinetic energy, i.e., sticking, for 
all v < v p . 

3.4 Combined coefficient of restitution 

The results from previous subsections can now be combined in Eq. (U6bb to compute 
the coefficient of restitution as a function of impact velocity for the irreversible elasto- 
plastic contact model presented in Fig. 4: 


6 — £o^n^i — 


( 1 ) 

Ej 

( 2 ) 

£/ 


for Vi < v p 
for Vi > v p 


(25) 


with v p from Eq. ( IT9t and £„ = 1 or < 0 for reversible and irreversible non-contact 
forces, respectively. Note that, without loss of generality, also other shapes of non- 
contact, possibly long-range interactions can be used here to compute e a and £,-, how¬ 
ever, going into these details goes beyond the scope of this paper, which only covers 
the most simple, linear non-contact force. 


4 Dimensionless parameters 

In order to define the dimensionless parameters of the problem, we first introduce the 
relevant energy scales, before we use their ratios further on: 


1 , 

Intial kinetic energy : £) = —m r vf , 


Attractive non — contact potential energy : E a = 


1 f 2 

_ Ja_ 

2 ki 


(26a) 


1 2 

Potential energy stored at <5£ ax : E p = -k\8^ ax , (26b) 


(26c) 


The first two dimensionless parameters are simply given by ratios of material parame¬ 
ters, while last two (independent) are scaled energies: 


Plasticity : rj = 


kp k\ 

k\ 


k 

Plastic (contact) adhesivity : B = — , 

k\ 


Non — contact adhesivity : a = 


f 2 

J a 


k\m r vj 


2 ’ 


Dimensionless (inverse) impact velocity : if/ = */ — = ^ max \l— 

V Ei Vi V m r 


(27a) 

(27b) 

(27c) 

(27d) 
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from which one can derive the dependent abbreviations: 


Scaled maximal deformation: 


Dimensionless impact velocity: 


X = 


8 P 

u ma 


V 
¥ v 


1 -fa¬ 


ct 


c = - = 


Vp v /t// 2 -2at// 


(28a) 

(28b) 


Using Eqs. (127ab . (I27bb and (128al) . ko can be rewritten in non-dimensional form 


Hx) 

k\ 

so that Eq. d22l> becomes: 


1 + riX, if Z < 1 
l + q, if % > 1 


(29) 


c n — 


(3t) 2 Z 2 


,l + t7Z (! + J7Z)(l+i3 + t7Z) 

and, similarly, Eq. (l24l > in non-dimensional form reads: 


y/-% 2 -2 ayx 


(30) 


e 


( 2 )_ 



P 7 ! 2 


(l + T7)(l+j3 + t7) 




(31) 


where the abbreviation y/X = 1 + a 2 + aj was used. 

To validate our analytical results, we confront our theoretical predictions with the 
results of two-particle DEM simulations in Fig. 5, which shows e plotted against the 
dimensionless impact velocity £, for elasto-plastic adhesive spheres with different non- 
contact adhesion strength f a (and thus a). The lines are the analytical solutions for the 
coefficient of restitution, see Eqs. (l30l > and (OH . and the symbols are simulations, with 
perfect agreement, validating our theoretical predictions. 

For low velocity, the coefficient of restitution e is zero, i.e., the particles stick to 
each other. This behavior is in qualitative agreement with previous experimental and 
numerical results li2Tl[54l[56 :l. With increasing impact velocity, e begins to increase 
and then decreases again, displaying a second sticking regime (for the parameters used 
here). For even higher impact velocity, v > v p (and thus £ > 1), another increase is 
observed, which will be explained in more detail in the next section. 

Besides the onset of the plastic-limit regime at £ = 1, we observe three further 
velocities ^, a \ and , for the end of the first sticking regime, i.e. 0 < £ < £c“\ 
and the second sticking regime, i.e. <C<d c) . While is constant, the critical 
velocity, required to separate the particles increases, whereas ^!' b \ required to 
enter the high-velocity sticking regime, decreases with the non-contact adhesion f a . 

The further study of these critical velocities and the comparison to existing litera¬ 
ture (theories and experiments) goes beyond the scope of the present study, since there 
are just too many possibilities for materials and particle sizes. We only refer to one 
example, where the end of the low-velocity sticking regime was predicted as a non¬ 
linear function of the surface energy/adhesion l26l and hope that our proposal of using 
dimensonless numbers will in future facilitate calibration of contact models with both 
theories and experiments. 
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Figure 5: Restitution coefficient e plotted as a function of the dimensionless impact 
velocity £, see Eq. (I28bb . for elasto-plastic adhesive spheres with the irreversible non- 
contact branch. Solid and dashed lines correspond to the analytical expressions in Eqs. 
(f30b and d3TT >. respectively, and the squares and circles are results of DEM simulations 
for values of f a as given in the legend in terms of a = 0.042 and 0.42 (for impact 
velocity v,- = 0.01 ms -1 ). Simulation parameters used here are k\ = 1() 2 Nm *, k p = 
5 x 10 2 Nm _1 , ( r/ = 4), k c = 10 2 Nm _1 , (/) = 1), with = 0.1, for particles with 
radius 1.1 x 10~ 3 m, density 2000kg/m 3 , and thus mass m r = 5.6 x 10~ 6 kg. 
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Figure 6: Restitution coefficient plotted as a function of the scaled initial velocity % 
for a collision without viscous and non-contact forces (a = 0). The solid red line 
corresponds to the analytical expression in Eq. ( l32l >. the dashed blue line to Eq. ( |33| >. 
the thin black line represents the low velocity approximation, and the circles are DEM 
simulation data. The material parameters are as in Fig. 5, i.e. r) = 4 and j8 = 1. 

5 Results using the mesoscopic contact model only (f a = 

0 ) 


Having understood the results for the contact model with finite non-contact force f a , we 
will restrict our analytical study to the special case of negligible non-contact adhesive 
forces f a = 0, in the following, which corresponds to the range of moderate to large 
impact velocity or weak f a , i.e. aCl. 

For this special case, the dimensionless parameters reduce to a = 0, and % = T = 
C = ¥-• The expressions for the coefficients of restitution presented in Eqs. (130b and 
(OTI reduce to 


and 


en\r],l5,x < 1 ) 

en\ri,p,x > 1 ) 




Pv 2 x 2 

(l + f7X)(l+j3 + T?z) 


n\ll,P,X = l ) 2 


(32) 

(33) 


5.1 Qualitative Description 

In Fig. 6, the analytical prediction for the coefficient of restitution, from Eqs. ( l32t 
and ( l33t . is compared to the numerical integration of the contact model, for different 
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scaled initial velocities %. We confirm the validity of the theoretical prediction for the 
coefficient of restitution in the whole range. 

For very small % < 10 2 , the approximation 1 ss 1 — ^ predicts the data very 
well. With increasing initial relative velocity, dissipation increases non-linearly with 
the initial kinetic energy, leading to a convex decrease of e„ * (due to the log-scale plot). 
The coefficient of restitution e ib becomes zero when a critical scaled initial velocity 
X, b] (see Eq. ( 1341 below) is reached. At this point, the amount of dissipated energy 
becomes equal to the initial kinetic energy, leading to sticking of particles. The coef¬ 
ficient of restitution remains zero until a second critical scaled initial velocity Xc is 
reached, i.e. sticking is observed for Xc^ ^ X < Xc' ■ Finally, for x > if 1 , the dis¬ 
sipated energy remains constant (the elastic limit branch is reached), while the initial 
kinetic energy increases. As a result, the kinetic energy after collision increases and 
so does the coefficient of restitution e n . Existence of sticking at such high velocities 
was recently reported by Kothe et al. (99], who studied the outcome of collisions be¬ 
tween sub-mm-sized dust agglomerates in micro-gravity can Note that an increase 
of e n for high velocity is a familiar observation in studies focused on the cold-spray 
technique CHU- Above a certain (critical) velocity the spray particles adhere to the 
substrate, and they do so for a range of impact velocities; increasing the impact ve¬ 
locity further leads to unsuccessful deposition, i.e. the particles will bounce from the 
substrate. The sticking and non-sticking phenomenology of such materials has been 
extensively studied experimentally and numerically in Refs. mm 

In Fig. 7, we compare the variation of the force with overlap in the various regimes 
of x as discussed above, but here for (j)f = 0.05. For very small %, the unloading slope 
k. 2 ~k\, (see Fig. 7(a) for a moderate X = 0.34), and the amount of dissipated energy is 
small, increasing with X- The kinetic energy after collision is almost equal to the initial 
kinetic energy, i.e. e„ ~ 1, see Fig. 6. In Figs. 7(b) and 7(c), the force-overlap variation 
is shown for sticking particles, for the cases Xc’' < X < 1 and 1 < X < Xc ' • respectively 
(more details will be given in the following subsection). Finally, in Fig. 7(d), the case 

(c) 

X > Xc is displayed, for which the initial kinetic energy is larger than the dissipation, 
resulting in separation of the particles. The corresponding energy variation is described 
in detail in annendixlBl 

5.2 Sticking regime limits and overlaps 

In this section we focus on the range of Xc’ 1 < Z < Xc'' ^ where the particles stick to 
each other (implying that /3 is large enough /3 > j3* with minimal j3* for sticking) and 
calculation of the critical values Xc^ and Xc' ■ When X = Xc” all of the initial kinetic 
energy of the particles is just dissipated during the collision. Hence the particles stick 
and <?!' ^ (rj, /3, x} b> ) = 0, which leads to 1 + j3 + r\X — P il 2 X 2 = 0- Only the positive 

7 Note that £ 1 is the regime where the physics of the contact changes, dependent on the material and 

other considerations; modifications to the contact model could/should then be applied, however, this goes 
beyond the scope of this paper, where we use the elastic limit branch or the generalized fully plastic model 
without it. Beyond the limits of the model, at such large deformations, the particles cannot be assumed to be 
spherical anymore and neither are contacts isolated from each other. 
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Figure 7: Contact force during one collision, plotted against the overlap for different 
scaled initial velocities x = 0.34, 0.69, 1.1, and 1.37, respectively. The three straight 
lines represent the plastic branch, with slope k\, the adhesive branch, with slope —k c , 
and the limit branch with slope k p , for k\ = 10 2 Nm“ 1 , k p = 5 x 10 2 Nm " 1 , k c = 10 2 
Nm _1 , i.e. r] = 4 and J3 = 1, and <j)f = 0.05. 
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solution is physically possible, as particles with negative initial relative velocity cannot 
collide, so that 


yiP) = J_ 

Xe 2)3/7 


1 + ^ 1 + 40 ( 1 + 0 ) 


(34) 


is the lower limit of the sticking regime. For larger x > Xc b \ the dissipation is strong 
enough to consume all the initial kinetic energy, hence the particles loose their kinetic 
energy at a positive, finite overlap 8 C , see Fig. 7(b). The contact deforms along the path 
0 — > <5 max -+ do —> <S mm —> 8 C . Thereafter, in the absence of other sources of dissipation, 
particles keep oscillating along the same slope k 2 . In order to compute 8 C , we use the 
energy balance relations in Eqs. < [20t . and conservation of energy along () mm — > 8 C , 
similar to Eq. (I20eb . 

- \k c { 8^ n - S 2 } = 0 , (35a) 

with vanishing velocity v c = 0 at overlap 8 C . Using the definitions around Eqs. (|2( )b 
and re-writing in terms of k c and 5 max leads to 


K-8 + t 1 - 




k c (k 2 -ki) 2 


k 2 k 2 (k 2 +k c ) 


5^ = 0 


(35b) 


and thus to the sticking overlap in regime (1), for 




ib) 


Vp < Vi < Vp\ 


8c l) _ g max / ( k 2 -k\) 2 k\ 
8U ~ 8U\ k 2 {k 2 + k c ) k 2 k c - 


(35c) 


In terms of dimensionless parameters, as defined earlier, this can be written as: 


5c (1) / 11 2 X 2 1 2 .(i) 

<S£ ax V ( l + r >X){l+P + rix) j3(l + ?7 x) v /^ £ ’" 


(36) 


where 1 denotes the result from Eq. (l32t with positive argument under the root. 

For larger initial relative velocities, X > 1, the coefficient of restitution is given by 
Eq. (l33b . so that the upper limit of the sticking regime Xc^ > 1 can be computed by 
setting e^'(r).fi.x'( C) ) =0. Again, only the positive solution is physically meaningful, 
so that 


yi C ) - 

X c — 



Pv 2 

(l + rj)(l+)3 + rj) 


(37) 


is the maximum value of x for which particles stick to each other. For X + X,' 1 particles 
deform along the path 0 —> 5m ax -+ 5 max — > 5o — > 5 rn in — > 8 C and then keep oscillating 
on the branch with stiffness k 2 , with 8 C being one of the extrema of the oscillation, 
see Fig. 7(c). Similar to the considerations above, we compute the sticking overlap in 

(c) 

regime (2), for v p < v; < Xc v p : in dimensionless parameters: 


4 (2) I V 2 . V X 2 _ X 42) 

<5£ ax ]/ (i + n)(i+P+v) P(1 + V) (3 yp e " 


(38) 
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Figure 8: Kinetic energy-free contact overlap <5 C plotted as a function of the scaled 
initial velocity x = v,/v |; ; the increasing branch corresponds to x < 1, while the de¬ 
creasing branch corresponds to x > 1- The dots are simulations for rj = 4 and J3 = 1, 
as in Fig. 7, which yields <S max /^max = (1 /3) l/2 in Eq. ( l39l >. 

where e„ 2) denotes the result from Eq. (l33l > with positive argument under the root. 

In Fig. 8, the scaled sticking overlap 8 C /S ^ax is plotted for different x, showing 
perfect agreement of the analytical expressions in Eqs. (l36t and < 138b . with the numeri¬ 
cal solution for a pair-collision. In the sticking regime, the stopping overlap increases 
with x, and reaches a maximum at x = 1, 

S C {X = !) = / /3r/ 2 — (1 + ?7 + J3) 

<5max y + J ?)(l + V + P) 

which depends on the the adhesivity /3 and the plasticity rj only. For x > 1, dissipation 
gets weaker, relatively to the increasing initial kinetic energy, and 5^ 2 ^/5£ ax decreases 

(c) 

until it reaches 0 for X = Xc ■ 

5.3 Contacts for different adhesivity p 

In the previous subsections, we studied the dependence of the coefficient of restitution 
e„ on the scaled initial velocity x for fixed adhesivity J3, whereas here the dependence 
of e n on ft is analyzed. 

A special adhesivity j3* can be calculated such that e n = 0 for x = 1, which is the 
case of maximum dissipation, and leads to sticking only at exactly x = 1, he- there is 
no sticking for /3 < j3*. Using Eq. (l32l) . we get 

1 +/3* + ?7 — P*rj 2 = 0 , (40a) 
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Figure 9: Coefficient of restitution <?„ plotted against the scaled initial velocity %. Cir¬ 
cles with different colors correspond to different adhesivity p (red for p < p* , green 
for /3 = J3* and blue for p > P*) for % < 1 , while magenta, black and cyan squares 
correspond to the respective values of p for x > 1 ■ Other parameters used are k\ = 10 2 , 
k p = 5 x 10 2 , and different k c (all in units of Nm" 1 ), i.e. rj = 4 and P/p* = 1/3, 1, 
and 3, with p* = 1 /3. The dashed red line represents the solution with the tuned fully 
plastic model with a new (j)f = 0.5 and newly calculated k p , see Appendix lAl. 


so that 

P* = ——r ■ (40b) 

V 1 

In Fig. 9, we plot the coefficient of restitution e as a function of the scaled initial 
velocity x for different values of adhesivity p. For p < /3*, in Fig. 9, the coefficient of 
restitution e„ decreases with increasing x < 1, reaches its positive minimum at x = 1, 
and increases for x > 1 ■ In this range, the particles (after collision) always have a non¬ 
zero relative separation velocity v/. When p = p*, e n follows a similar trend, becomes 
zero at x = 1, and increases with increasing scaled initial velocity for x > 1- This is 
the minimum value of adhesivity for which e„ can become zero and particles start to 
stick to each other. For p = p*, the sticking regime upper and lower limits coincide, 
Xc^ = xf ] = 1. If P > P*, e n decreases and becomes zero at x = Xc '* < 1, it remains 
zero until x = Xc > 1, arid increases with increasing relative initial velocity thereafter. 
Hence, the range of velocities for which sticking happens is determined by the material 
properties of the particles. Indeed Zhou el al. E) presented similar conclusions about 
the deposition efficiency in cold spray. Simulations with viscous forces change the 
value of P* and are not shown here, see Appendix [B] 
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6 Conclusions and outlook 


Various classes of contact models for non-linear elastic, adhesive and visco-elasto- 
plastic particles are reviewed. Instead of focusing on the well understood models for 
perfect spheres of homogeneous (visco) elastic or elasto-plastic materials, a special 
mesoscopic adhesive (visco) elasto-plastic contact model is considered, aimed at de¬ 
scribing the macroscale behavior of assemblies of realistic fine particles (different from 
perfectly homogeneous spheres). An analytical solution for the coefficient of restitu¬ 
tion of pair-contacts is given as reference, for validation, and to understand the role of 
the contact model parameters. 

Mesoscopic Contact Model The contact model by Luding fl5l . including short- 
ranged (non-contact) interactions, is critically discussed and compared to alternative 
approaches in subsec tion ll. 21 The model introduced in section[2]is simple (piece-wise 
linear), yet it catches the important features of particle interactions that affect the bulk 
behavior of a granular assembly, i.e. non-linear elasticity, plasticity and contact adhe¬ 
sion. It is mesoscopic in spirit, i.e. it does not resolve all the details of every single 
contact, but is designed to represent an ensemble of particles with many contacts in a 
bulk system. One goal of this study is to present this rich, flexible and multi-purpose 
granular matter meso-model, which can be calibrated to realistically model ensembles 
of large numbers of particles HlOOl . The analytical solution for the contact dissipation 
is given for contact and non-contact forces both active, but viscosity inactive, in sec¬ 
tion [3] A sensible set of dimensionless parameters is defined in section [4j before the 
influence of the model parameters on the overall impact behavior is discussed in detail, 
focusing on the irreversible, adhesive, elasto-plastic part of the model, in section [5] 

Analysis of the coefficient of restitution When the dependence of the coefficient 
of restitution, e, on the relative velocity between particles is analyzed, two sticking 
regimes, e = 0, show up, as related to different sources of dissipation: 

(i) As previously reported in the literature (see e.g. Refs. Il2ll|52ll55ll56l[62l ) the 
particles stick to each other at very low impact velocity. This can happen due to irre¬ 
versible short-range non-contact interaction, as e.g. liquid bridges, or due to van der 
Waals type force for dry adhesive particles. The threshold velocity, below which the 
particles stick depends on the magnitude of the non-contact adhesive force f a , while 
for elasto-plastic adhesive particles on both non-contact adhesive force and plasticity, 
which together control this low-velocity sticking. 

(ii) With increasing velocity, e increases and then decreases until the second stick¬ 
ing regime is reached, which is strongly influenced by the plastic/adhesive (and vis¬ 
cous) dissipation mechanisms in the hysteretic contact meso-model. At small impact 
velocity, all details of the model are of importance, while at higher velocities, for a suf¬ 
ficiently low value of the jump-in force f a , the contribution of the non-contact forces 
can be neglected. The theoretical results are derived in a closed analytical form, and 
phrased completely in terms of dimensionless parameters (plasticity, adhesivity and 
initial velocity). The ranges of impact velocities for the second sticking regime are 
predicted and discussed in detail. 
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(iii) For still increasing relative velocity, beyond the second sticking region, e starts 
increasing again. This regime involves a change of the physical behavior of the sys¬ 
tem as expected, e.g. for non-homogeneous materials with micro-structure and non-flat 
contacts, or materials with an elastic core, e.g. asphalt (stone with bitumen layer). Even 
though this elastic limit behavior is a feature of the model, completely plastic behavior 
can be reproduced by the model too, just by tuning two input parameters k p and 0y, as 
shown in appendix[A] This way, the low velocity collision dynamics is kept unaffected, 
but the elastic limit regime is reached only at higher impact velocities, or can be com¬ 
pletely removed. This modification provides the high velocity sticking regime for all 
high velocities, as expected for ideally plastic materials. On the other hand, the exis¬ 
tence of a high velocity rebound, as predicted the model with elastic limit regime, has 
been observed experimentally and numerically in cold spray (2}{6l and can be expected 
for elastic core with a thin plastic shell. 

Additional dissipative mechanisms For sticking situations, on the un-/re-loading 
branch, the particles oscillate around their equilibrium position until their kinetic en¬ 
ergy is dissipated, since realistic contacts are dissipative in nature. Since viscosity 
hinders analytical solutions, it was not considered before, but a few simulation results 
with viscosity are presented in Appendix[B] With viscosity, both un- and re-loading are 
not elastic anymore, resembling a damped oscillation and eventually leading to a static 
contact at finite overlap. 

Application to multi-particle situations The application of the present meso-model 
to many-particle systems (bulk behavior) is the final long-term goal, see e.g. Ref. l28l 
|29][35), as examples, where the non-contact forces were disregarded. An interesting 
question that remains unanswered concerns a suitable analogy to the coefficient of 
restitution (as defined for pair collisions) relevant in the case of bulk systems, where 
particles can be permanently in contact with each other over long periods of time, and 
where impacts are not the dominant mode of interaction, but rather long lasting contacts 
with slow loading-unloading cycles prevail. 

One specific example for the latter situation of slow loading-unloading of bulk 
material is given in Appendix[F] showing qualitatively similar behavior as encompassed 
in the contact meso-model, but on a much larger length-scale than the contact model 
itself, highlighting the dominant role of the particle structure and the (non-flat) contact 
area with related plastic (irreversible) re-arrangements |26| . 

Outlook The interest of widely different communities (viz. particle technology, gran¬ 
ular physics, interstellar dust, asphalt, or cold-spray) in the dependence of the coeffi¬ 
cient of restitution (or deposition/impact behavior) on the impact velocity is consid¬ 
erable. We hope our study helps to connect these widely different communities by 
providing an overview and, in particular, a flexible, multi-purpose contact model, valid 
and useful for many practically relevant situations. 

The contact meso-model has to be calibrated for different materials, while our refer¬ 
ence analytical results allow to verify the model implementation. With this, the model 
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can be used to predict bulk material behavior and to be validated by comparison with 
experiments. 
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A Tuning of parameters to increase the plastic range 

We assume that the reference dimensionless plasticity depth be <j>f , which is, e.g. cal¬ 
culated based on the maximal volume fraction related arguments of a multi-particle as¬ 
sembly, and k p be the reference limit stiffness. We propose a new (j)f > (j)f, which rep¬ 
resents the new (larger) dimensionless plasticity depth (arbitrary choice or calculated 
based on another volume fraction consideration) and a new value of k p ; the choice is 
such that the tuned model resembles exactly, consistently, the reference for So < c/i?©/, 
with reduced radius a 12 , and becomes plastic for a^tj)/ < Sq < .At 5a = ai20/, 
Eq.® reads 

kp = h + {k;~h)8P i J8^, (A.l) 

since all parameters except and k p remain unchanged. Using Eqs. dA. 1 b and (fill we 
arrive at 

(k p — k\) 2 _ (k p ~k\) 2 
k p <t>f k p <j)f 

which gives the new limit stiffness 

k p = k x +AB/2 + ^(AB/2) 2 + kxAB, 

where A = (k p —k\) 2 /k p and B = 0y /0y. 

Using Eq. (1A.3I >. we can calculate values of the new limit plastic stiffness k p for 
any given <j)f , such that the collision dynamics for lower plastic deformation <5o < 8 P 0 
is intact, while the range of plastic deformation is enhanced, depending on the chosen 
<j)f ><t> f . 


(A.2) 

(A.3) 


B Effect of Viscosity 

Since real physical systems also can have additional dissipation modes that are, e.g. 
viscous in nature, we study the behavior of collisions with viscous damping present 
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Figure 10: Force-displacement law for elasto-plastic, adhesive contacts superimposed 
on the irreversible contact force law. The black solid line represents the force law for 
reference input parameters (j)r and k p , while the dashed red line represents the same for 

I / 

a new chosen (j)/ and newly calculated k p resembling a wider plastic regime of the 
particle deformation. 
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(yo > O) and compare it with the non-viscous case (yo = 0). Note that any non-linear 
viscous damping force can be added to the contact laws introduced previously, how¬ 
ever, for the sake of simplicity we restrict ourselves to the simplest linear viscous law as 
given as second term in Eq. (0, since it can be important to choose the correct viscous 
damping term for each force law to get the physically correct behavior P6ll50ll74]|l0ll . 
In Fig. 11, we plot the contact force against the overlap, and the overlap against time, 
during collisions for a constant value of x = 1 and different /3, for yo = 5 x 10 3 kg s' 1 . 

When ft < p*, see Fig. 11(a) and Fig. 11(b), the contact ends when the adhesive 
force — k c 8 goes back to zero, for both cases, with and without viscosity. This is since 
the viscosity is relatively small and does not contribute enough to the total dissipation 
to make the particles stick for the parameters used here. 

For the critical adhesivity p = ft*, reported in Fig. 11(c), without viscosity, the 
overlap between the particles goes down to exactly zero at the end of the collision, 
with all kinetic energy dissipated. For yo > 0, dissipation brings this marginal case into 
the sticking regime and the particles stay in contact at 8 > 0. This can be seen clearly 
in Fig. 11(d), where the particles undergo a damped oscillatory motion due to the small 
residual velocity created on the re-loading branch. 

For larger values p > /3*, the overlap at which kinetic energy is lost completely (on 
the k c branch) is finite, for both yo = 0 and yo > 0, see Fig. 11(e). In both cases, the 
particles stick and remain in contact. Without viscosity, the particles keep oscillating 
along the slope kj, while with viscosity the oscillation is damped and kinetic energy 
vanishes. During loading and unloading the apparent slope changes with time due to 
the additional viscous force that leads to the dissipation of energy, as evident from the 
ellipsoidal converging spiral. Waiting long enough, for some oscillation cycles, the 
particles stick to each other with a finite overlap and zero relative kinetic energy. The 
difference is also visible in Fig. 11(f), where for yo = 0 the particles keep oscillating 
with constant amplitude, whereas, for yo > 0, the particles undergo a damped oscilla¬ 
tory motion, until the velocity becomes 0 at 5 > 0. The time evolution of the overlap in 
Fig. 11(f) resembles that of the displacement evolution in Ref. j 102] , where the authors 
studied sticking of particles in Saturn’s rings. @ 


C Asymptotic Solutions 

In this subsection, we focus on the case j < 1, and study the asymptotic behavior of 
the coefficient of restitution as function of the impact velocity. 

For the sake of simplicity, let us start with an elasto-plastic system without adhe¬ 
sion, i.e. k c = 0, in Eq. (l32l > such that 


e£\ri,p=0,x<l ) = yYT^’ (Cla) 

8 In general, one could add a viscous law that is proportional to &2 — &i or to a power of overlap <5, such 
that the jump-in viscous force in Fig. 11(e) at the beginning of the contact is not there, however, we do not 
go into this detail. 
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Figure 11: (a), (c), (e) Contact forces plotted against overlap and (b), (d), (f) time 

evolution of 5/c>max for pair collisions with parameters k\ = 10 2 , k p = 5 x 10 2 and 
different k c = 10, 33.33, and 100, (units Nm -1 ), i.e. with 77 = 4, (5 < /3*, /3 = /3* and 
f5 > j5*, for the same situations as shown in Fig. 9. The red and blue lines represent the 
data in the presence and absence of viscosity respectively, where yo = 5 x 10 3 , (unit 
Nm^'sec). 






































Figure 12: The coefficient of restitution is plotted against the scaled initial velocity 
X in log-log-scale for /3 = 0 and three values of Tj = 5, 50, and 500, with the other 
parameters as in Fig. 6 . Red, green and blue circles denote, respectively, the solution 
of Eq. Q, while the solid lines represent the approximation for high scaled impact 
velocity and large plasticity t] 1 . 


inserting the definitions of 17 , /3 and v p , 


ei X \p = 0 ,v< v p ) 


1 


N 


1+ 


kp k\ Vi 



(C.lb) 


using Eq. (ITB . where we defined S = f 1 ’/ 1 and assuming ca 0 = \/ —-, we get 

M^max V 171 

e { n\p = 0,V < Vp) = Aj—^7- (C.lc) 

V 1 + 4 

Eq. (OB is exactly the same as Eq. (5) in fZD- For non-cohesive particles, and in the 
range v < v p we get exactly the same solution as Walton and Braun 1711 . 

Further to study the asymptotic solution 


ei 1 ) (t7,/3=0,x<l) = y'^^«(t 7 ^) (C.2) 

with the approximation valid for i]j> 1. Since the scaled velocity is moderate, % < 
1 , the condition requires a large plasticity, i.e., a strong difference between the limit 
stiffness and the plastic loading stiffness, q 1 (or k p k\). In Fig. 12, we plot the 
coefficient of restitution against the scaled initial velocity % for three different values 
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Figure 13: Log-log plot of the coefficient of restitution against the scaled initial velocity 
X for four different values of p = 0.01, 0.1, and 1.0, with r] =50. Red, green and blue 
circles denote the respective solutions of the general equation, Eq. ( l32l >. solid black line 
represents power law e„ ~ V l/4 , while magenta line denotes e„ ~ v -1 / 2 . 


of t] = kp/k i, together with the power law prediction of Eq. (1C.21) . We observe, that 
for the smallest rj (red circle and line), the approximation is far from the data, while 
for higher rj , the approximation works well even for rather small velocities X ~ 0 . 1 . 

Next, when studying the elasto-plastic adhesive contact model, p > 0 and p <C 1, 
again, we restrict ourselves to values of 77 such that asymptotic condition r/x ^ 1 is 
satisfied. Hence, Eq. (l32l > can be approximated as 


e { n\ri,P,x < 1 ) 



(C.3) 


as long as r\x ^ P > 0 and ^ > p holds. 

In Fig. 13, we plot the coefficient of restitution against the scaled initial velocity x 
for different values of P and superimpose the approximation, Eq. (1C.31) . For small p 
and large j, one observes good agreement between the full solution and the approxi¬ 
mation. Differently, for the highest values of P the approximation is not valid. Due to 
the adhesive force, for large x> with increasing p, the deviation from the X 1 2 power 
law becomes increasingly stronger, leading to the sticking regime, as discussed in the 
previous subsections. On the other hand, for smaller velocities, one observes a con¬ 
siderably smaller power-law, resembling the well-known j -1 / 4 power law for plastic 
contacts, as indicated by the dashed line in Fig. 13. 
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Figure 14: Semi-log plot of the coefficient of restitution as function of the scaled initial 
velocity using different interpolation rules for ki, for pair collisions with T) = 50 
and P = 0. The symbols denote the solutions of the general equation, Eq. (1C.2I > with 
linear interpolation (red circles) or square root interpolation (blue circles), as given in 
Eq. dED. The red and blue solid lines represent the approximations for high impact 
velocity e„ ~ X~ 1 ^ 2 an d e « ~ Z ^ 4 - 

D Dependence on interpolation 

The choice of the interpolation rule for the unloading stiffness £2 in Eq. ( 129} is empiri¬ 
cal. Therefore, for 5 max /5max < 1, a different choice could be: 

*2(4nax)=* 1 (l+TJv/z). (D.l) 

Inserting Eq. (ID. II ) into Eq. (OTT ) leads to a different expression for the normal coeffi¬ 
cient of restitution ei ] \ which for high values of T) and for small /), reduces to 

e n x -v/t?(z) _1/4 ■ (D.2) 

A similar power law prediction for moderate velocities has been previously obtained 
by Thornton et al. in Ref. ED, using a non-linear Hertzian loading and unloading. Fig. 
14 shows the agreement between the power law approximation x 1/<4 and Eq. (|2H with 
the alternative interpolation rule (ID. 1 b . for moderate velocities. The choice of different 
interpolation laws for ki shows the flexibility of the model and requires input from 
experiments to become more realistic. The convexity of linear interpolation for zero 
cohesion is very similar to that of low P in Fig. 9. 
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E Energy Picture 


This appendix shows the energies of two particles during contact, where the difference 
between the different branches of the contact model, namely irreversible/unstable or 
reversible/elastic, will be highlighted. 

In Fig. 15, the time-evolution of kinetic and potential energy is shown; the graphs 
can be viewed in parallel to Figs. 7(a) and 7(b). In Fig. 15(a), we plot the kinetic 
and potential energy of the particles against time for low initial velocity % < X^ ! ■ 
corresponding to Fig. 7(a), for which dissipation is so weak that particles do not stick. 
The kinetic energy decreases from its initial value and is converted to potential energy 
(the conversion is complete at 5 max ). Thereafter, the potential energy drops due to the 
change between the loading and unloading slope from k\ to ki. The potential energy 
decreases to zero (at the force-free overlap <5o), where it is converted to (less) kinetic 
energy. Then the kinetic energy decreases further due to the acting adhesive force. 
At <5 mm the increasing potential energy drops to a negative value due to the change in 
unloading slope from k 2 to the adhesive (instable) slope —k c . From there it increases 
from this minimum, negative value to zero, for 5 = 0. From here the kinetic energy 
remains constant and the potential energy stays at zero, since the particles are separated. 

In Fig. 15(b), we plot the time evolution of kinetic and potential energy that the par¬ 
ticles would have if un-/re-loading would take place at that moment, along the branch 
with slope k 2 , namely the available (elastic) potential energy. This energy increases 
from zero at t = 0, and reaches a maximum when the kinetic energy becomes zero 
(note that it is not equal to the initial kinetic energy due to the plastic change of slope 
of k 2 -) Thereafter, the available potential energy decreases to zero at the force-free 
overlap 5 (l . For further unloading, the available potential energy first increases and 
then drops rapidly on the unstable branch with slope — k c . The change in sign of the 
unloading slope, from ki to —k c , is reflected in the kink in the curve at 5 m j n . Note, that 
comparing Figs. |15(a)| and [f5(b)[ the available potential energy always stays positive, 
while the total, plastic “potential” energy drops to negative values after the kink at 5 m j n . 

Figs. 15(c) and 15(d) show the time evolution of kinetic and potential energy (total 
and available, respectively) for an initial velocity Xc^ < X < X^ i n the sticking regime, 
see Fig. 7(b). In Fig. 15(c), a similar trend as that of Fig. 15(a) is observed until the 
potential energy becomes negative at 5 mm . The difference to the case of smaller impact 
velocity is that at this point, the kinetic energy is less than the magnitude of the negative 
potential energy and hence first reaches zero, i.e., the particles stick. At this point, the 
(plastic) potential energy increases and jumps to a positive value indicating the change 
in sign of the unloading slope from —k c to /jo. Finally, it oscillates between this positive 
value at 5 C , exchanging energy with the kinetic degree of freedom. When the available 
potential energy is plotted in Fig. 15(d), a similar trend as that of Fig. 15(b) is observed 
up to the kink at 5 m j n . Here, the two energies have comparable values when they reach 
5 m ; n and the kinetic energy decreases to zero with a non-zero available potential energy, 
which causes the contact to re- and un-load along £ 2 - 
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(a) 



(c) 



(b) 



(d) 


Figure 15: (a), (c) Kinetic and (irreversible, plastic, “potential”) energy of the particles, 
and (b), (d) kinetic and available (elastic) potential energy (for re-loading) of the parti¬ 
cles, plotted against time for pair collisions with k\ = 10 2 Nm -1 . k p = 5 x 10 2 Nm~ 1 , 
and k c = 10 2 Nm _1 , i.e. 77 = 4 and j 8 = 1. The initial velocity % is % = 0.34 (a,b) and 
X = 0.69 (c,d), in the regimes defined in the inset of each plot. 
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F Cyclic agglomerate compression and tension tests 


Goal of this appendix is to show the unloading and re-loading behavior of an agglomer¬ 
ate, i.e. its effective, mesoscopic force-displacement relation, which clearly is different 
from the contact force law applied at the primary particle contacts. We will report 
incomplete detachment and partly/weaker elastic response for re-loading after various 
different compressive and tensile loading amplitudes. 

The system considered here is an agglomerate (cubic) of size Lq =0.115, made of 
N = 1728 primary particles of diameter do = 0.01 (with some variation in size to avoid 
monodisperse artefacts), just as in Ref. ESI. The cubic sample was first compressed 
(pressure-sintered) with a dimensionless wall stress dop s /k p = 0.02 to form a stable, 
rather dense agglomerate or “tablet”. The stress is first released to a value 2.10 5 , i.e. 
Pr/Ps = 10 for all walls. Then various uni-axial, unconfined tension/compression 
tests are carried out applying either further tension or compression starting from the 
released state of the sample El. The simulation parameters are same as in Ref. ESI 
(table 2), except for the cohesion that is set here to a rather small intensity, k c /k p = 0.2, 
rolling and sliding friction coefficients that are double as large, /I, = ji„ = 0.2, and 
viscous damping of those degrees of freedom, y r /y = y„/'y = 0.1, which also is larger 
than that of the reference situation. 

The force-displacement curves for the tests at different amplitudes are shown in 
Figs. 16 and 18 for tension and compression tests respectively. All simulations in 
Figs. 16 and 18 start from the same configuration, i.e. the released state mentioned 
above and is indicated by the black circle at point (0,0). These plots represent the 
mesoscopic contact model of agglomerates consisting of multiple primary particles 
and their geometrical surface configurations and change in shape during the tests. 

Fig. 16 shows the force-displacement curve for an unconfined uniaxial tension test. 
The black arrow shows the unloading/tension path, and finally arrows with different 
colors show the re-loading paths for different deformation amplitudes, as given in the 
inset. Each of the tests, when it reaches the original strain at zero, is then repeated for 
three more cycles. Note that repeated cyclic loading remains on the same branch with 
positive slope, displaying the elastic nature of the contact, while it is not completely, 
perfectly detached. The contact surface is changing plastically by restructuring of the 
primary particles and surely is not flat, see Fig. 17, as one would expect for ideal, ho¬ 
mogeneous, plastic materials. For the largest amplitude, the behavior is not perfectly 
elastic anymore, since the first plastic effects show up. For deformations as large as 0.2 
of the primary particle diameter, do, before re-loading (arrow with positive slope on 
the red curve) has mostly, but not completely lost mechanical contact. The complete 
detachment of the assembly happens for much higher amplitude, than what is expected 
from a two-particle interaction. Note that the contact model of the primary particles 
is behaving elasto-plastically (<pr = 0.05) on the scale of only 0.05iio; the reversible, 
elastic un-/re-loading is thus not due to the primary particle contact model, since it 
stretches to four times Of do and even higher displacements. Finally, in order to con¬ 
firm that this is not an effect of viscosity, qualitatively, the thick lines are simulations 
performed four times slower than those with thin lines. 

In Fig. 17, a few snapshots during the tensile deformation are presented. The first 
snapshot corresponds to the undisturbed sample, while the others are increasing tensile 
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Figure 16: Dimensionless force-displacement curve for an unconfined uni-axial ten¬ 
sion test (negative horizontal axis), with the various different deformation amplitudes 
D x given in the inset. The downward arrow indicates the direction of first tensile un¬ 
loading, while the upwards-right arrows indicate the change of force during re-loading. 
Except for the red curve, all these branches are reversible, for repeated un-/re-loading. 



Figure 17: Snapshots of the tablet-sample during (large) tensile deformations for D x = 
( L-L 0 )/d 0 = 0 (a), 0.81 (b), 1.8 (c), 3.1 (d), 4.7 (e), 7.4 (f), and 8.6 (g). The primary 
particles are colored according to their distance from the viewer (red, green, blue is 
increasing distance). 
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(L-L 0 ) / d 0 


Figure 18: Dimensionless force-displacement curve for the same sample as in Fig. 
16, but under compressive initial loading and un-/re-loading. The values in the inset 
indicate the maximal amplitudes D x . 


deformation amplitudes. Note that these deformations are much larger than in Fig. 16. 
The contact is completely lost only at the extreme, final deformation in Fig. 19(g). In 
Fig. 17, it is also visible that the contact surface has developed a roughness of the size 
of several primary particles; the first visible gap is opened at a total deformation of 
D x ~ do, and the contact is lost only at D x ~ 8c/o, when the last of the thin threads 
breaks. The elastic, irreversible tension branch, however, is strongly developed only 
for much smaller D x ~ do/5. 

Complementing the tension test above. Fig. 18 shows the behavior of the same sam¬ 
ple during compression cycles. The values given in the inset indicate the amplitude of 
un-/re-loading. The smallest amplitudes remain elastic throughout, while plastic defor¬ 
mation kicks in forD* >0.1 (see the red curve). However, the unloading and re-loading 
take place on the same branch, i.e. a new elastic branch (e.g. for D x = 0.2). For even 
larger amplitudes, e.g. the yellow curve with D x = 0.3, the continuous damage/plastic 
destruction of the sample (by considerable irreversible re-arrangement during each cy¬ 
cle). Again, thick lines indicate simulations four times slower, which shows a small 
quantitative difference, but qualitative agreement even for the largest amplitude/rate. 
The snapshots in Fig. 19 show the continuous plastic deformation of the sample at 
large strains. 
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